clear;
load('gr_data_all_20230413_v2.mat');
close all
clc
a = gr_data_all;
T = unique(a(:,1));

MGcolor = [235, 147, 52]/255;

G = [];
Gerr = [];
for k=1:6
    g = reshape(a(:,k+1),[length(a)/length(T) length(T)]);
    gmean = nanmean(g,1)';
    gstd = nanstd(g,1)';

    G = [G gmean];
    Gerr = [Gerr gstd];

end


%% Plot data
T1plot = 16;
T2plot = 50;

clrs = distinguishable_colors(6);
clrs = cbrewer('qual', 'Paired',6);
strainnames = {'MG1655', 'BW25113', 'CS109', 'BL21', 'dnaK', 'uup'};

ms = 6;
f1 = figure;

toplot = [1];

for k=toplot
    g_plot = G(:,k);
    gerr_plot = Gerr(:,k);
%     errorbar(T,g_plot,gerr_plot,'ko', 'markersize', ms, 'markerfacecolor', MGcolor);
    errorbar(T,g_plot,gerr_plot,'ko', 'markersize', ms, 'markerfacecolor', [112, 150, 212]/255);
    hold on;
    
end

MG_growthrate = g_plot;
MG_growthrate_err = gerr_plot;

set(gcf, "Position", [0 0 400 300]);
set(gca, 'FontSize', 20);
xlabel('Temperature (°C)');
ylabel('Growth rate (1/h)');
box off;
legend(strainnames{toplot})
xlim([T1plot T2plot])

%% Define fit variables
Tfit = T(4:11);
gfit = MG_growthrate(4:11);
gfit_err = MG_growthrate_err(4:11);

%% Perform linear fit and plot fit
figure(f1);
hold on;
flm_linear = fitlm(Tfit, gfit, 'weights', gfit_err);
plot(Tfit, flm_linear.feval(Tfit), 'k','linewidth', 1)
rss = sum(abs(gfit - flm_linear.feval(Tfit)).^2)
flm_linear

%% Perform Arrhenius fit
Arrhenius_fit = @(b,x) b(1).*exp(-b(2)./(x+273.15));
b0 = [1E12 8000];

flm_Arrhenius = fitnlm(Tfit, gfit, Arrhenius_fit, b0, 'Weights', gfit_err);
hold on;
plot(Tfit, flm_Arrhenius.feval(Tfit), 'r', 'linewidth', 1)
rss = sum(abs(gfit - flm_Arrhenius.feval(Tfit)).^2)
flm_Arrhenius

%% Perform Eyring fit
% Eyring_fit = @(b,x) b(1)./(x + 273.15) .*exp(-b(2)./(x+273.15));
% b0 = [1E12 8000];
% 
% flm_Eyring = fitnlm(Tfit, gfit, Eyring_fit, b0, 'Weights', gfit_err);
% hold on;
% plot(Tfit, flm_Eyring.feval(Tfit), 'g', 'linewidth', 1)
% flm
